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Abstract —We consider the problem of estimating the distance, 
or range, between two locations by measuring the phase of a 
sinusoidal signal transmitted between the locations. This method 
is only capable of unambiguously measuring range within an 
interval of length equal to the wavelength of the signal. To 
address this problem signals of multiple different wavelengths 
can be transmitted. The range can then be measured within 
an interval of length equal to the least common multiple of 
these wavelengths. Estimation of the range requires solution of 
a problem from computational number theory called the closest 
lattice point problem. Algorithms to solve this problem require 
a basis for this lattice. Constructing a basis is non-trivial and 
an explicit construction has only been given in the case that the 
wavelengths can be scaled to pairwise relatively prime integers. In 
this paper we present an explicit construction of a basis without 
this assumption on the wavelengths. This is important because the 
accuracy of the range estimator depends upon the wavelengths. 
Simulations indicate that significant improvement in accuracy 
can be achieved by using wavelengths that cannot be scaled to 
pairwise relatively prime integers. 

Index Terms —Range estimation, phase unwrapping, closest 
lattice point 


I. Introduction 


Range (or distance) estimation is an important component 
of modern technologies such as electronic surveying 01 
and global positioning 0 1 . Common methods of range 
estimation are based upon received signal strength [5}, 6], 
time of flight (or time of arrival) 0 80, and phase of ar¬ 
rival 0 EL m This paper focuses on the phase of arrival 
method which provides the most accurate range estimates in 
many applications. Phase of arrival has become the technique 
of choice in modern high precision surveying and global 
positioning 


m 


A difficulty with phase of arrival is that only the principal 
component of the phase can be observed. This limits the 
range that can be unambiguously estimated. One approach to 
address this problem is to utilise signals of multiple different 
wavelengths and observe the phase at each. Range estimators 
from such observations have been studied by numerous au¬ 
thors 01 [30. Least squares/maximum likelihood and 
maximum a posteriori (MAP) estimators of range have been 
studied by Teunissen 141], Hassibi and Boyd [12], and more 
recently Li et. al. 0. A key realisation is that least squares 
and MAP estimators can be computed by solving a problem 
from computational number theory known as the closest lattice 
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point problem 113, 14]. Teunissen 0 appears to have been the 
first to have realised this connection. 

Efficient general purpose algorithms for computing a closest 
lattice point require a basis for the lattice. Constructing a 
basis for the least squares estimator of range is non-trivial. 
Based upon the work of Teunissen 0 , and under some 
assumptions about the distribution of phase errors, Hassibi and 
Boyd II1211 construct of a basis for the MAP estimator. Their 
construction does not apply for the least squares estimator!] 
This is problematic because the MAP estimator requires 
sufficiently accurate prior knowledge of the range, whereas 
the least squares estimator is accurate without this knowledge. 
An explicit basis construction for the least squares estimator 
was recently given by Li et. al. 10 under the assumption 
that the wavelengths can be scaled to pairwise relatively prime 
integers. In this paper, we remove the need for this assumption 
and give an explicit construction in the general case. This is 
important because the accuracy of the range estimator depends 
upon the wavelengths. Simulations show that a more accurate 
range estimator can be obtained using wavelengths that are 
suitable for our basis, but are not suitable for the basis of 
Li et. al. ifToh . 

The paper is organised as follows. Section |TT] presents the 
system model and defines the least squares range estimator. 
Section m introduces some required properties of lattices. 
Section m shows how the least squares range estimator is 
given by computing a closest point in a lattice. An explicit 
basis construction for these lattices is described. Simulation 
results are discussed in Section [V] and the paper is concluded 
by suggesting some directions for future research. 


II. Least squares estimation of range 
Suppose that a transmitter sends a signal x(t) = e 2 ^Ut+v) 
with phase tj> and frequency / in Hertz. The signal is assumed 
to propagate by line of sight to a receiver resulting in the 
signal 

y(t) = ax(t — ro/c) + w(t) = ae 2n ^ t+e ' ) + w(t ) 

where r o is the distance (or range) in meters between receiver 
and transmitter, c is the speed at which the signal propagates 
in meters per second, a > 0 is the real valued amplitude of 
the received signal, w{t) represents noise, 0 = (ft — vq/\ is the 
phase of the received signal, and A = c// is the wavelength. 

1 The least squares estimator is also the maximum likelihood estimator under 
the assumptions made by Hassibi and Boyd El- The matrix G in El is rank 
deficient in the least squares and weighted least squares cases and so G is 
not a valid lattice basis. In particular, observe that the determinant of G E3 
p. 2948] goes to zero as the a priori assumed variance cr^ goes to infinity. 
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The receiver is assumed to be synchronised by which it is 
meant that the phase f and frequency / are known to the 
receiver. 

Our aim is to estimate ro from the signal y(t). To do this 
we first calculate an estimate 9 of the principal component 
of the phase 9. In optical ranging applications 8 might be 
given by an interferometer. In sonar or radio frequency ranging 
applications 9 might be obtained from the complex argument 
of the demodulated signal y(f)e _27r M Whatever the method of 
phase estimation, the range ro is related to the phase estimate 
9 by the phase difference 

Y=(J>-e) = (ro/\ + $} ( 1 ) 

where $ represents phase noise and ( x) = x — \_x + 4J where 
[xj denotes the greatest integer less than or equal to x. For 
all integers k, 

Y = ( r o/A + $) = ((ro + k\)/\ + $}, (2) 


is called an n-dimensional lattice. The matrix B is called a 
basis or generator for A. The basis of a lattice is not unique. 
If U is an n x n matrix with integer elements and determinant 
det U = ±1 then U is called a unimodular matrix and B and 
BU are both bases for A. The set of integers 7L n is called 
the integer lattice with the n x n identity matrix I as a basis. 
Given a lattice A its dual lattice, denoted A*, contains those 
points that have integral inner product with all points from A, 
that is, 

A* = {x ; x'y € Z for all y € A}. 


The following proposition follows as a special case of Propo¬ 
sition 1.3.4 and Corollary 1.3.5 of 1181. 


Proposition 1. Let v G Z n , let H be the n — 1 dimensional 
subspace orthogonal to v, and let 


vv 

Q = I—— =I ~ 

■\r' 


vv' 

IMF 


and so, the range is identifiable only if ro is assumed to lie 
in an interval of length A. A natural choice is the interval 
[0, A). This poses a problem if the range ro is larger than the 
wavelength A. To alleviate this, a common approach is to trans¬ 
mit multiple signals x„(t) = e 27r ^ nt+ ^ for n = 1 
each with a different frequency f n . Now N phase estimates 
9\,... ,9n are computed along with phase differences 

Y n = {(/)- 9 n ) = (r 0 /A n + $„) n = 1,... ,N (3) 


where A n = c/f n is the wavelength of the nth signal 
and represent phase noise. Given Yi,..., Y N , a 

pragmatic estimator of the range ro is a minimiser of the least 
squares objective function 


N 

LS(r) = Y J {Y n ~ r/K } 2 . 

n =1 


(4) 


This least squares estimator is also the maximum likelihood 
estimator under the assumption that the phase noise variables 
$i,..., <I> v are independent and identically wrapped normally 
distributed with zero mean [15, p. 501 lfl6l p. 76lll7[ p. 47], 
The objective function LS is periodic with period equal to 
the smallest positive real number P such that P/X n G Z for 
all n = 1,..., N, that is, P = lcm(Ai,...,A n) is the least 
common multiple of the wavelengths. The range is identifiable 
if we assume ro to lie in an interval of length P. A natural 
choice is the interval [0, P) and we correspondingly define the 
least squares estimator of the range ro as 


f = arg min LS(r). (5) 

r£[0,P) 

If A n /Afc is irrational for some n and k then the period P 
does not exist and the objective function LS is not periodic. 
In this paper we assume this is not the case and that a finite 
period P does exist. 


III. Lattice theory 

Let B be the mxn matrix with linearly independent column 
vectors bi,..., b„ from m-dimensional Euclidean space R m 
with m > n. The set of vectors 


be the n x n orthogonal projection matrix onto H. The set 
of vectors Z" fl H is an n — 1 dimensional lattice with dual 
lattice (Z” fl H)* = {Qz ; z e Z“}. 


Given a lattice A in M m and a vector y £ R m , a problem of 
interest is to find a lattice point x £ A such that the squared 
Euclidean norm ||y—x|| 2 = Y^iLi(Vi~ x i) 2 ' s minimised. This 
is called the closest lattice point problem (or closest vector 
problem) and a solution is galled a closest lattice point (or 
simply closest point) to y 0. 

The closest lattice point problem is known to be NP- 
hard |T^, 2(J. Nevertheless, algorithms exist that can compute 
a closest lattice point in reasonable time if the dimension is 
small (less that about 60) 114, 2~iT- FUl. These algorithms have 
gone by the name “sphere decoder” in the communications 
engineering and signal processing literature. Although the 
problem is NP-hard in general, fast algorithms are known for 
specific highly regular lattices [25l [2611 . For the purpose of 


range estimation the dimension of the lattice will be N — 1 
where N is the number of frequencies transmitted. The number 
of frequencies is usually small (less than 10) and, in this case, 
general purpose algorithms for computing a closest lattice 
point are fast [14]. 


IV. Range estimation and the closest lattice 

POINT PROBLEM 

In this section we show how the least squares range esti¬ 
mator f from (01 can be efficiently computed by computing a 
closest point in a lattice of dimension iV — 1. The derivation is 
similar to those in 1271 - 10 . Our notation will be simplified by 
the change of variable r = P/3, where P is the least common 
multiple of the wavelengths. Put v n = P/A„ £ Z and define 
the function 

N 

P(/3) = LS(Pf3) = J2 CM - K) 2 ■ 

n= 1 

Because LS has period P it follows that F has period 1. 
If /3 minimises F then P/3 minimises LS and, because f £ 
[0, P), we have f = P(f3 — |_/3J). It is thus sufficient to find a 
minimiser /3 £ 1 of F. 


A = {Bu ; u G Z”} 
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Observe that (Y n - /3v n ) 2 = min zeZ (Y n - / 3v n 

so F may equivalently be written 


z) 2 and 


our construction to be simpler than that in 0 - The following 
proposition is required. 


N 

F((3) = min VVlji - fiv n - z n ) 2 . 

n= 1 

The integers Zi,... ,Zn are often called wrapping variables 
and are related to the number of whole wavelengths that 
occur over the range r 0 between transmitter and receiver. The 
minimiser /3 can be found by jointly minimising the function 

N 

Fi(P,zi,...,z n ) = Y n - fiv n - z n ) 2 

n= 1 

over the real number f3 and integers z\,... ,Zn- This minimi¬ 
sation problem can be solved by computing a closest point in 
a lattice. To see this, define column vectors 

y = (Y 1 ,...,Y N ) / &R n , 

z = (zi,.. ,,z N y g z N , 

V = v N y = (P/Xu ■ ■ •, P/Xn)' G Z n . 

Now 


F 1 (P,z 1 ,...,z n ) =F 1 {/3 ,z) = ||y-/3v-z|| 2 . 

The minimiser of Fi with respect to ji as a function of z is 


Proposition 2. Let U be an N x N unimodular matrix with 
first column given by v. A basis for the lattice A* is given by 
the projection of the last N — 1 columns of U orthogonally 
onto H. That is, Qu 2 ,.... Qu y is a basis for A* where 
Ui,..., u/v are the columns of U. 

Proof: Because U is unimodular it is a basis matrix for 
the integer lattice Z N . So, every lattice point z G Z N can be 
uniquely written as z = ciUi + - • -+cnUn where ci,..., cjv G 
Z. The lattice 

A* = {Qz ; z G Z N } 

= {Q(ciUi + • • • + cat u at) ; ci,..., cat G Z} 

= {c 2 Qu 2 + • • • + catQujv ; c 2 ,..., cat G Z} 

because Qui = Qv = 0 is the origin. It follows that 
Qu 2 , .... Quat form a basis for A*. ■ 

To find a basis for A* we require a matrix U as de¬ 
scribed by the previous proposition. Such a matrix is given 
by Li et. al. [10, Eq. (76)] under the assumption that the 
wavelengths can be scaled to pairwise relatively prime in¬ 
tegers. We do not require this assumption here. Because 
P = lcm(Ai,..., A at) it follows that the integers tq,..., Vn 
are jointly relatively prime, that is, gcd(rq,..., vn) = 1 . 
Define integers <q,..., g at by gw = vn and 


Substituting this into Fj gives 

F 2 (z) = minF 1 (/3, z) = F 1 (j3(z),z) = ||Qy- Qz || 2 

/3eR 

where Q = I — vv'/||v|| 2 is the orthogonal projection matrix 
onto the N — 1 dimensional subspace orthogonal to v. Denote 
this subspace by H. By Proposition Q] the set A = Z w n H is 
an N — 1 dimensional lattice with dual lattice A* = {Qz ; z G 
Z N }. We see that the problem of minimising F 2 (z) is precisely 
that of finding a closest point in the lattice A* to Qy G M' v . 
Suppose we find x G A* closest to Qy and a corresponding 
z G Z N such that x = Qz. Then z minimises F 2 and /3(z) 
minimises F. The least squares range estimator in the interval 
[0, P) is then 

r = P(j3( z)-L/3(z)J). (6) 

It remains to provide a method to compute a closest point 
x G A* and a corresponding z G Z ;Y . In order to use known 
general purpose algorithms we must first provide a basis for 
the lattice A* 0 . The projection matrix Q is not a basis 
because it is not full rank. As noted by Li et. al. [IQ], a 
modification of the Lenstra-Lenstra-Lovas algorithm due to 
Pohst [30] can be used to compute a basis given Q. However, it 
is preferable to have an explicit construction and Li et. al. IllOh 
give a construction under the assumption that the wavelengths 
Ai,..., Aat can be scaled to relatively prime integers, that is, 
there exists cGi such that gcd(cAfc, cA„) = 1 for all k no 
We now remove the need for this assumption and construct a 
basis in the general case. As a secondary benefit, we believe 

-The assumption that the wayelengths are pairwise relatively prime is made 
implicitly in equation (75) in G3- 


g k = gcd(v fe , ...,v N )= gcd(v k ,gk+i), k = 1,..., N - 1 


and observe that gk+i/gk and Vk/g k are relatively prime 
integers. For k = 1,..., N — 1, define the A by A matrix 
A k with m, nth element 


Vk/gk 

gk+i/gk 


Akmn — * 




b k 


m = n = k 
m = k + 1 , n = k 
m = k,n = k + 1 
m = n = k + 1 

otherwise 


where I mn = 1 if m = n and 0 otherwise. The integers a k 
and bk are chosen to satisfy 


bk 


Vk_ 

gk 


- a k 


fffc+i 

gk 


= 1 


(7) 


and can be computed by the extended Euclidean algorithm. 
The matrix A/, is equal to the identity matrix everywhere 
except at the 2 by 2 block of indices k < m < k + 1 and 
k < n < k + 1. The matrix A k is unimodular for each k 
because it has integer elements and because the determinant 
of the 2 by 2 matrix 


Vk/gk a k _ , v k g k +i _ , 

/ 7 — Ok CLk — J- 

gk+i/gk o k g k gk 

as a result of ®. A matrix U satisfying the requirements of 
Proposition [2] is now given by the product 


N -1 

U = Ak = An-i x Ajv -2 x ■ • • x Ap 

fc=l 
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Fig. 1. Comparison of the least squares range estimator with wavelengths 
A (top) and C (bottom) suitable for the basis in o and B (top) and D 
(bottom) suitable only for the basis described in this paper. Sets B and D 
result in smaller mean square error when the noise variance cr 2 is greater than 
approximately 1.2 x 10 —4 and 7 x 10“ 5 respectively. 
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That U is unimodular follows immediately from the unimod¬ 
ularity of Ai,..., Ajv-i- It remains to show that the first 
column of U is equal to v. Let vr,..., vjy-i be column 
vectors of length N defined as 

v fc = (v 1 ,...,v k ,g k+1 ,0,...,0)', k = 1,..., AT-2 
vjv_i = (ui,.. .,v N -!,g N y = v. 


One can readily check that v/, :+ i = A/,. + i v/. for all k = 
1,..., N — 1. The first column of the matrix Ai is Vi and so, 
by induction, the first column of the product Ylk=i A/,. is 
for all K = 1,..., N — 1. It follows that the first column of 
U is vjv-i = v as required. 

Let U 2 be the N by N — 1 matrix formed by removing 
the first column from U, that is, U 2 = ( 112 ,..., ujv). By 
Proposition[2]a basis for A* is given by projecting the columns 
of U 2 orthogonally to v, that is, a basis matrix for A* is the 
N by N — 1 matrix B = QU 2 . Given B a general purpose 
algorithm [14] can be used to compute w £ T, N ~ l such that 
x = Bw is a closest lattice point in A* to Qy £ M. N . Now 


x = Bw = QU 2 W = Qz 

and so z = U 2 W £ Z N . The least squares range estimator f 
is then given by (O. 


V. Simulation Results 

We present the results of Monte-Carlo simulations with 
the least squares range estimator. Simulations with N = 4 
and N = 5 wavelengths are performed. For each case we 
consider two different set of wavelengths. The first set is 
suitable for the basis of Li et. al. 01 and was used in the 
simulations in [10]. The second set is suitable only for our 
basis. In each simulation the true range r ( ] = 20 and the phase 
noise variables $ 1 ,..., <&jv are wrapped normally distributed, 
that is, <l> n = (X n ) where X t ..... X N are independent 
and normally distributed with zero mean and variance cr 2 . In 
this case, the least squares estimator is also the maximum 
likelihood estimator. Figure |T| shows the sample mean square 


error for cr 2 in the range 10 -5 to 10 -2 and 10 7 Monte-Carlo 
trials used for each value of cr 2 . 

For N = 4 the two sets of wavelengths are 

A _ f 9 q r ij \ td _ f 210 210 210 210 ~] 

Ai JD — \ yg , 61 , 41 , 31 /• 

For both sets the wavelengths are contained in the interval 
[2, 7] and P = 210 = lcm(A) = lcm(B) so that the identifi¬ 
able range is the same. The wavelengths A are relatively prime 
integers and are suitable for the basis of Li et. al. Hi and are 
used in the simulations in 11 Oil. The wavelengths B are not 
suitable for the basis of llOll because they can not be scaled 
to pairwise relatively prime integers. To see this, observe that 
the smallest positive number by which we can multiply the 
elements of B to obtain integers is c = 61 ^q 49 - Multiplying 
the elements of B by c we obtain the set 

c x B = {77531,100409,149389,197579} 


and these elements are not pairwise relatively prime because, 
for example, gcd(77531,100409) = 1271. Figure |T| shows 
the results of simulations with both sets A and B. When the 
noise variance a 2 is small wavelengths A result in slightly 
reduced sample mean square error as compared with B. 
As a 2 increases the sample mean square error exhibits a 
‘threshold’ effect and increases suddenly. The threshold occurs 
at cr 2 ~ 1.2 x 10~ 4 for wavelengths A and cr 2 ss 3 x 10~ 4 
for wavelength B. Wavelengths B are more accurate than A 
when a 2 is greater than approximately 1.2 x 10 -4 . 

For N = 5 the two sets of wavelengths are 

_ /9 ^ ^ 7 111 D -T 2310 2310 2310 2310 2310 \ 

^ 1 5 - Ly l 877 ’ 523 ’ 277 5 221 ’ 211 J * 

For both sets all wavelengths are contained in the interval 

[2,11] and P = 2310 = lcm(C') = lcm(U) so that 
the maximum identifiable range is the same. The basis of 
Li et. al. J3 can be used for wavelengths C but not for 
D. The wavelengths C were used in the simulations in E3i. 
Figure Q] shows the result of Monte-Carlo simulations with 
these wavelengths. Wavelengths C result in slightly smaller 
sampler mean square error than D when cr 2 is small, but 
dramatically more error for a 2 above the threshold occurring 
at (J 2 « 7 x 10" 5 . 

The sets B and D have been selected based on a heuristic 
optimisation criterion. The properties of this criterion are not 
yet fully understood and will be the subject of a future paper. 


VI. Conclusion 

We have considered least squares/maximum likelihood es¬ 
timation of range from observation of phase at multiple 
wavelengths. The estimator can be computed by finding a 
closest point in a lattice. This requires a basis for the lattice. 
Bases have previously been constructed under the assumption 
that the wavelengths can be scaled to relatively prime integers. 
In this paper, we gave a construction in the general case and 
indicated by simulation that this can dramatically improve 
range estimates. An open problem is how to select wavelengths 
to maximise the accuracy of the least squares estimator. We 
will study this problem in future research. 
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